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Summary 

A long-standing, though ill-understood problem in rocket dynamics, rocket response to ran- 
dom, altitude-dependent nozzle side-loads, is investigated. Side loads arise during low altitude 
flight due to random, asymmetric, shock-induced separation of in-nozzle boundary layers. In 
this paper, stochastic evolution of the in-nozzle boundary layer separation line, an essential 
feature underlying side load generation, is connected to random, altitude-dependent rotational 
and translational rocket response via a set of simple analytical models. Separation line motion, 
extant on a fast boundary layer time scale, is modeled as an Ornstein-Uhlenbeck process. Pitch 
and yaw responses, taking place on a long, rocket dynamics time scale, are shown to likewise 
evolve as OU processes. Stochastic, altitude-dependent rocket translational motion follows from 
linear, asymptotic versions of the full nonlinear equations of motion; the model is valid in the 
practical limit where random pitch, yaw, and roll rates all remain small. Computed altitude- 
dependent rotational and translational velocity and displacement statistics are compared against 
those obtained using recently reported high fidelity simulations [Srivastava, Tkacik, and Keanini, 
J. Applied Phys., 108, 044911 (2010)]; in every case, reasonable agreement is observed. As an 
important prelude, evidence indicating the physical consistency of the model introduced in the 
above article is first presented: it is shown that the study's separation line model allows direct 
derivation of experimentally observed side load amplitude and direction densities. Finally, it 
is found that the analytical models proposed in this paper allow straightforward identification 
of practical approaches for: i) reducing pitch/yaw response to side loads, and ii) enhancing 

pitch/yaw damping once side loads cease. 
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1. Introduction 

Although the recorded history of powered rocket flight spans a millenium, dating to eleventh century 
China [Ij, and experimentally-based rocket development traces two hundred years to eighteenth 
century England and while an enormous scientific and engineering literature attaches to the 
dynamics, design, and control of rockets, missiles, and spacecraft, numerous interesting, practically 
important questions remain. One of the most basic and well-studied concerns prediction of rocket 
trajectories during ascent. While the long time success of large rocket and missile programs suggests 
that this essential question has been solved, in fact, a variety of difficult-to-predict features introduce 
significant uncertainty. These include: i) altitude-, attitude-, and speed-dependent aerodynamic 
forces [Ij, ii) random loads produced by location- and altitude-dependent wind and atmospheric 
turbulence (3j|6], iii) rocket design and construction imperfections [7j, iv) fuel sloshing (liquid fuel 
rockets) [8|[9], v) slag formation (solid fuel rockets) [2], and vi) random impacts with air-borne 



animals and debris 10 



Random, nozzle side loads present a further, singularly ill-understood feature complicating as- 



cent prediction. As depicted in Fig. 1, and as described, e.g., in 11,12 , nozzle side loads appear 
in over-expanded nozzles during low- altitude flight, when high ambient pressure forces external 
air upstream into the nozzle. Under these conditions, the inflow overcomes the low pressure, low 
inertia, near-wall nozzle out-flow. At a locus of points, the nominal, instantaneous boundary layer 
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separation line, inflow inertia decays sufficiently that outflow inertia can turn the inflow back on 
itself; the reversed inflow forms a virtual corner and, in turn, an oblique, circumferential shock. See 
Fig. 2. 



Oblique 
Shock 




Figure 1: Shock- induced boundary layer separation in overexpanded supersonic nozzle flow. The 
process typically occurs during low altitude flight when ambient pressure is high enough to force 
atmospheric air into the nozzle. The incoming air flows upstream along the low-inertia, near-wall 
region until downstream-directed boundary layer inertia turns it, forming a virtual compression 
corner. An oblique shock thus forms, and the combined action of shock-induced pressure rise and 
inertial pressurization produced by the inflow forces the down-flow boundary layer to separate. 
Pressures, mach numbers, and velocities are denoted, respectively, by P, M, and U and V. Axial 
positions where the boundary layer starts to thicken (i denotes incipient), and where it separates 
are denoted, respectively, as Xi and Xs] the nominal shock-boundary layer interaction zone is shown 
as Lg. Since the separation line position, Xg, and downstream conditions vary with the altitude- 
dependent ambient pressure. Pa = Pa{H{t)) [12] , all variables shown likewise vary with H{t). 



Random side loads appear due to the pressure jump, AP, extant across the asymmetric, stochas- 
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tically evolving separation line. Viewing the instantaneous separation line as the superposition of 
the line's slowly moving mean axial position, Xs{t), and it's rapid random fluctuation about Xs{t), 
one observes that since Xs{t) and AP are both altitude-dependent, random side loads, and resulting 
torques and translational and rotational rocket responses, are all likewise altitude-dependent. 




Oblique 
Shock 



Separation 
Line 



Figure 2: Schematic of stochastic boundary layer separation line and associated, rippled, azimuthal 
oblique shock. The mean separation line position relative to the nozzle throat, Xg, varies with rocket 
altitude, H{t); the corresponding nozzle radius is i? = R{H{t)). The instantaneous separation line 
position relative to Xs{t) is shown as s{(j),t). The separation line lies on the nozzle wall and, in 
a nominally symmetric nozzle, the shock forms an azimuthally independent, average angle which 
varies with Xs{t). 



Supersonic flow separation and side load phenomena in nozzles, studied under fixed exter- 



nal pressure, i.e., effectively fixed altitude conditions, has attracted significant attention 113-22 



Ostlund and Muhammad-Klingmann l23| review much of this work. Numerical and experimental 
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studies have been reported, for example, by Frey and Hagemann (15,24 , Onofri and Nasuti 25 , 



Pekkari |26|, Schwane and Xia |27|, and Shimizu et al. f28'. Fundamental studies of shock-boundary 



layer interactions [29-34 in which (typically oblique) shocks are produced by a variety of shock 
generators, have, in turn, provided essential insight into the fluid dynamic features underlying 
shock-induced separation. 

Nozzle side loads can be perilous, damaging nozzles and/or attached equipment, or inducing 
catastrophic excitations within structural components; the failure of the Japanese H-II in 1999, for 
example, has been attributed to nozzle side loads |35|. In large engines, side load magnitudes can 
be extreme, on the order of 250,000 pounds, for example, in Saturn V engines (7). Thus, minimizing 
and designing to accommodate side loading represents an essential rocket design task. 

While side loads have long complicated ascent prediction, until recently [IT], their effects on 
rocket dynamics have remained completely unknown. The work in [11| focused on the stochastic 
ascent of sounding-rocket-scale rockets subject to altitude-dependent random side loads and torques. 
A high fidelity numerical model was developed which incorporated the following features: 

a) a variable mass, six-degree-of-freedom, nonlinear rocket dynamics model |36|, 



b) empirical models of the altitude-dependent ambient atmosphere [5] and mach number-dependent 
drag coefficient [37], 



c) a model of separation line motion which incorporated experimentally observed statistical 



properties of shock-separated flat plate turbulent boundary layers [31 32 34 , and 
d) an ad hoc, though physically consistent model of side load statistics. 
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The model allowed Monte Carlo simulation of random rocket ascents, as well as estimates of 
altitude-dependent, ensemble- averaged translational and rotational velocities and displacements. 



1.1. Objectives 

Through the remainder of the paper, we will often refer collectively to the set of numerical models 
described in |11| as Model /; the present set of analytical models will likewise be referred to as 
Model II. 

A) The first, and most crucial objective centers on establishing the physical consistency, i.e., 
the congruence with physical principles and available experimental data, of both Models 
I and II. Due to the difficulties associated with detecting and isolating dynamic side loads 
within nozzles subject to simultaneous random aerodynamic and structural loading, obtaining 
experimental data on the effects of side loading on ascent will likely remain problematic. 
Indeed, no such data presently exists. 

Thus, three key features motivate pursuit of this objective: 

i) substantial experimental challenges will likely continue to limit direct measurements. 



ii) theoretical insight is presently limited to the numerical models comprising Model I 11 
and crucially, 

iii) demonstration of physical consistency establishes a reasonable foundation for further work 
on this long-standing problem. 

The first objective is pursued in two steps. 
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(a) We first show that the separation line model introduced in allows derivation of 
both the assumed side load model in [ll], as well as experimentally observed side load 
amplitude and direction densities [l8]420|. See section 2. 



Since Model I [11] incorporates a high fidelity rocket dynamics model and reasonable 
models of rocket aerodynamics and altitude-dependent ambient atmospheric conditions, 
and since the above demonstration ties the most uncertain portions of Model I, i.e., the 
separation line and side load models, to experimental observations [18| - |20| , we argue that 
Model I [llj thus represents a physically consistent description of rocket ascent under 
altitude-dependent side loads. 

(b) Presuming the physical consistency of Model I ^llj, we then argue that since Model II 
predicts altitude-dependent rotational and translational rocket-response statistics that 
are, in every instance examined, consistent with those estimated via Model I (see section 
7), Model II is likewise consistent. 

B) It is found that the simple models comprising Model II provide a straightforward framework 
for interpreting observed random rocket responses to side loading, as well as for identifying 
practical approaches for designing against side loads. See section 7. The second objective 
thus centers on highlighting and exploiting this theoretical framework. 

1.2. Overview 

The following models and interconnecting elements tie stochastic separation line motion to rocket 
response, and thus comprise the paper's essential frame: 
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i) Local, short-time-scale separation line dynamics are modeled as an Ornstein-Ulilenbeck pro- 
cess. (See section 3.) Here, as depicted in Fig. 2, local refers to axial separation line motion in 
the vicinity of any given in-nozzle azimuthal angle, (j), while the short time scale corresponds 
to the correlation time for local axial separation line displacements. We will sometimes refer 
to the latter as the boundary layer time scale. 

ii) Since side load evolution takes place on this short time scale, while rocket response to side 
loading occurs on a much longer time scale, it becomes necessary to derive a coarse grained 
time correlation function for local separation line displacement appropriate to the longer 
scale. This procedure in turn leads to the crucially important (long-time-scale) side load time 
correlation function. See section 4. 

iii) Given the latter, and focusing first on rocket rotational dynamics, we show in section 5 that 
the stochastic evolution of pitch and yaw likewise correspond to Ornstein-Uhlenbeck processes. 
Once this key result is obtained, then theoretical expressions for altitude-dependent means and 
variances of the pitch and yaw rate and pitch and yaw displacement follow. Additionally, and 
as becomes apparent when interpreting results in section 7, parametric relationships obtained 
for associated (effective) damping and diffusion coefficients provide practical guidance for 
minimizing the effects of side loads on rotational dynamics. 

iv) An asymptotic model of translational rocket motion, appropriate in the practical limit where 
pitch, yaw, and roll all undergo small random variations about zero, follows in section 6. 
Again, the simplified linear model, derived from the general nonlinear model in allows 
analytical determination of altitude-dependent means and variances for the rocket's lateral 
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velocity and displacement components. 

Once a complete model is thus established, altitude-dependent variances of pitch and yaw rate, 
pitch and yaw displacement, lateral velocities, and lateral displacements are computed and com- 
pared against ensemble average estimates obtained via the model in [TI] (section 7) . 

2. Physical consistency of Model I 

As detailed in flT\, and in response to the decaying altitude-dependent ambient pressure, the mean 
position of boundary layer separation line, Xs{t), travels down the nozzle axis toward the nozzle exit, 
with motion taking place on a relatively slow time scale, = Axa/Vji, where Axa is the charac- 
teristic incremental altitude over which ambient pressure varies and Vr is the characteristic rocket 
speed. Superposed on this slow motion is a fast, random, azimuthally homogeneous stochastic 
motion. Following flT], the joint probability density, ps, associated with the instantaneous random 
separation line shape is given by 

/ N TT ^ r Si -|- Sn ~l~ 'S'? ~l~ • ■ • ~l~ 'S^r^ /-,\ 

Psisi, S2, . . . , sn) = [[pi = (^2Tras)^/^ ''''^ ^ 2^1 ^ ^ ^ 

where, as shown in Fig. 3, sj is the random axial displacement of the separation line at azimuthal 
angle (/>/, and is the (assumed) constant variance of local separation line displacements. 

Constituent displacements in the set of N displacements are assumed independent, and based 
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on experimental observations [sT 32,34 , gaussian. Thus, each pj is given by 



1 



In moving to a continuous description of the separation Hne, (TTl assumes that 
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<sicP,t)s{c^',t)>=a^J{cP -(!)') (3) 

Considering next the side load, we express the instantaneous force vector produced by asym- 
metric boundary layer separation, Fs(i), as a sum of radial and axial components 

Fs(t) = Fr(t) + Fx(t) (4) 



In 11 , the following ad hoc side load model was assumed: 

A) Fgy and Fsz are independent, gaussian random variables, 

B) < Fsy >= and < Fsz >= 0, 

C) (^{Fsy- < Fsy > )') = ((F,,- < Fsz > )') = 

where, assuming ergodicity, < • > denotes either an ensemble or time average, and where the 
separation line model above is used to calculate the force variance u^. The side load components 
Fsy and Fsz are expressed with respect to rocket-fixed coordinates; see Fig. 4 in section 5 below. 

In order to demonstrate the physical consistency of Model I, we first show that the assumed 
properties. A) -C), can be derived from the simple separation line model in [IT]. We then follow 111] 
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and show that the model in A) - C) leads to experimentally observed side load amplitude and 
direction densities. 




Figure 3: Model I separation line model. The mean separation line position, Xs{H{t)), moves 
down the nozzle axis, on the slow time scale associated with vertical rocket motion. By contrast, 
axial separation line motion about Xs{H{t)), at any angular position, (pi, is random, and takes 
place on a much shorter time scale. Rapid axial motion, in addition, is confined to the nominal 
shock-boundary layer interaction zone, again denoted by Lg. Pressures upstream and downstream 
of the instantaneous separation line, Pi = Pi{H{t)) and P2 = P2{H{t)), respectively, are assumed 
to be spatially uniform within Lg. Adapted from 



12 



Considering assumption B) first, we calculate: 



n2-K 

{Fsy{H{t)))^ = R[H{t)][Pi{H{t))-Pa{H{t))\ / sin<P{s{^,t))d<i> (5) 

JO 

and 

/•27r 

{F,,{H{t)))^ = R[H{t)][P,{H{t))-Pa{H{t))\ / cos<p{s{(t>,t))dct> (6) 

Jo 

where R{H{t)) is the nozzle radius at the axial location of the mean separation line at time i, 
Pi{H{t)) — Pa{H {t)) = AP {H{t)) , is the pressure jump across the separation-inducing shock 
at altitude H{t), Pi and Pa are the wall pressures immediately upstream and downstream of the 
instantaneous separation line, and where we have approximated the downstream pressure as the 
instantaneous ambient pressure, Pa{H{t)) [3,8,15]. See [llj for further discussion. 

In order to evaluate these averages, express the k^^ realization of, e.g., Fgz, in discrete form as 

M 

S2,..., sm) = R[H{t)] [P^{H{t)) - Pa{H{t))] ^ s^''\<Pi) COS (t>iA^ (7) 

1=1 

where s^'^^cpi) is the associated separation line displacement at (/>/. Taking the ensemble average 
term by term, and noting that 

/oo poo poo 
/ .../ SlPs{si,S2, . . ■ ,SN)dSidS2- ■ ■ dSN = 
-oo J~oo J —oo 

then confirms B): 

{F,yiH))^ = (F,,(F)>^ = (8) 
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Turning next to A), since Fgy and Fgz represent sums of say, independent gaussian random 
variables, then both are gaussian. Computing 



< Fsy{H)Fs,iH) >= R\H)^P^{H)al [ J sin2</.# = 

Jo ^ 

where ([s]) has been used, Ukewise shows that Fsy{H) and Fsz{H) are independent. 

Finally, for C), since < F^y{H) >= R'^ (H) AP^iH)a'^ f^^ cos"^ (j)d(l), and 
< Fl{H) >= R^{H)AP'^{H)a1 " sin^ then 

< Fl{H) >=< Fl{H) > (9) 

2.1. Derivation of density functions for side load amplitude and direction 

We note two important experimental and numerical observations concerning the side load, Fg 
(within rigid, axisymmetric nozzles): 



a) the probability density of the random amplitude, A = |Fs|, is a Rayleigh distribution 18-20 
and 

b) the random instantaneous direction, (ps, of Fg is uniformly distributed over the periphery of 



the nozzle, or p^^{(j)s) = l/2vr, where p^p^ is the pdf of the side load direction 18-20 



As noted, both observations can be derived, starting from the simple statistical model of ran- 
dom side loads. A) - C), immediately above. Thus, given A and 0s, the instantaneous side load 
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components in body-fixed y and z directions (see Fig. 4) are given by 



Fsy = A cos ( 



F^z = ^sin( 



Following 11 , write Fsy and Fsz as Fsy = Y = Acos(f)s and Fsz = Z = Asm(j)s; thus, the joint 



probability density associated with Fsy and Fsz can be expressed as 



PYziY, Z) = Py{Y)pz{Z) = exp 



2^2 



(10) 



Following |llj, we restate pyz iii terms of A and (ps as, 



(11) 



where Pa<i>s 4's) is the joint pdf for the random amplitude and direction of mathbfFr, and where 
the jacobian determinant is given by 



\J\ 



dY dV 

dA d(t>s 

dZ dZ 

dA d<t>s 



A 



[12) 



Thus, 



PA<j,s{^,4>s 



A 



27rc72 



exp 



1 

2^ 



A 

exp 



'2cj2 



s)pa{A) (13) 
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where, 



P<^si<Ps) = ^ O<0, <27r (14) 



is the uniform probabihty density underlying the random direction (pg, and 



A / A'^ \ 

Pa{A) = -^expl-^j (15) 



is the Rayleigh distribution for the amphtude A. 

In summary, we have shown that the separation hue model in [11] : i) allows derivation of 
the assumed side load model. A) - C), in ^llj, and ii) provides a theoretical basis for explaining 



experimentally observed 18-20 side load amplitude and direction densities. 



3. Ornstein-Uhlenbeck model of separation line dynamics 

Theoretical determination of rocket response to side loads requires that the time correlation function 
for either side load component, {Fsa{t')Fsa{t)) , be first determined. This and the next section 
develops {Fsa{t')Fsa{t)) in two steps. First, we propose (and physically justify) that local separation 
line dynamics can be modeled as an Ornstein-Uhlenbeck process. Once this assumption is made, 
then the second step rests on a rigorous argument showing that on the relatively long rocket 
dynamics time scale, the boundary layer separation line shape, and importantly, associated side 
load components, are all delta correlated in time. See section 4. 
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We propose the following simple, explicit stochastic model of separation line dynamics: 



dsi{t) = -ksi{t) + y/WsdWit) (16) 

where Si (t) = s {(j)i,t) is the instantaneous separation line position at 4>i, k and Ds are damping 
and effective diffusion coefficients, and dW{t) is a differential Weiner process. This equation, 
describing an Ornstein-Uhlenbeck process, allows straightforward, physically consistent calculation 
of statistical properties associated with separation line motion and, more importantly, serves as the 
first link in a chain that connects short-time-scale random separation line motion to short-time-scale 
random side loads, and in turn, to long-time-scale stochastic rotational rocket dynamics. 

The form of this equation is chosen based on the following experimental features, observed in 
shock-separated flows near compression corners and blunt fins: 

a) Under statistically stationary conditions, the feet of separation-inducing shocks oscillate 
randomly, up- and downstream, over limited distances, about a fixed mean position; see. 



e.g., 31,32,34 . 



b) As observed in j32[|34| the distribution of shock foot positions within the shock-boundary 
layer interaction zone is approximately gaussian. 

c) The time correlation of shock foot positions, as indicated by wall pressure measurements 
within the shock-boundary layer interaction zone, decays rapidly for time intervals, At, larger 



than a short correlation time, Tg, a feature that can be inferred, for example, from 38 



Physically, the damping term captures the fact that the shock sits within a pressure-potential 
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energy well. Thus, downstream shock excursions incrementally decrease and increase, respectively, 
upstream and downstream shock face pressures; the resulting pressure imbalance forces the shock 
back upstream. A similar mechanism operates during upstream excursions. Introduction of a 
Weiner process models the combined random forcing produced by advection of turbulent boundary 
layer structures through the upstream side of the shock foot and pressure oscillations emanating 
from the downstream separated boundary layer and recirculation zone. 

We note that the proposed model is qualitatively consistent with Plotkin's model of boundary 



layer-driven shock motion near compression corners and blunt fins 38 . Plotkin's model, which 
captures low frequency spectra of wall pressure fluctuations within these flows, corresponds to a 
generalized Ornstein-Uhlenbeck process in which a deterministic linear damping term is superposed 
with a non-Markovian random forcing term. We use an ordinary OU process model, incorporating 
a Wiener process, since again, it is consistent with the above observations and more particularly, 
since it allows much simpler calculation of statistical properties. 

4. Derivation of the coarse grained side load time correlation func- 
tion 



Given the model in Eq. ( 16 ), our path shifts to obtaining a mathematically consistent description of 



resulting side load statistics. To accomplish this, we take advantage of the significant separation in 



time scales that exists between large-scale, low frequency random separation line motion 31 ,32 34 
and the relatively slow dynamics of the rocket. Thus, define the long time scale as i = t e~^, where 
e << 1 is a dimensionless scale factor and t is the (short) time scale associated with low frequency. 



18 



large scale separation line motion. The long time scale can be chosen to correspond to any of 
a number of dynamical features; the chosen scale determines e. Since we are interested in the 
dynamical response of the rocket, we choose i to be on the order of tr = LV^^, where L and Vr 
are the length and characteristic speed of the rocket. Thus, e = f^^T^^ , where / is a characteristic 



frequency from the low frequency band associated with large scale shock foot motion 31 ,32 . Based 
on the rocket parameters given in fTT], r^j is on the order of 2(10~^) s; thus, since is on the 
order of 10~^ s 



32 



34 



, e = 0(10-2) to 0(10-1^ 



Using well known expressions 39 for the mean and variance of an OU process, one can readily 
show that the process becomes stationary on time scales that are long relative to the characteristic 
period for large scale separation line motion, i.e., for t ^ k~^. Under stationary conditions, the 
variance becomes independent of time, var (si (t)) = Dg / {2k) , and the time correlation can be 
placed in the form: 

< Si (r) Si (0) >,= ^ exp(- |t| / Tc ) (17) 

where t = t2 — h, and Tc = < Si (t) Sj (0) > dt/ var (sj (t)) , is the correlation time. Under 
stationary conditions, Tc = and var (sj (t)) = DgTc / 2 = a^. 

Given the above time correlation, we can now show that on long time scales, t, < Si (r) Si (0) >«, 
approaches delta- function like character. First, define a normalized time correlation function as 



RAr) = ^^^^^^^^^ (18) 
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where Rs (r) dT=l. Next, define a rescaled version of Rg (r) as Rge (t) = e ^i?s (e ^r) , or 



i?..(r) = ^^if^^^ (19) 

Z(J „ T^£^ 



where r = e ^t. Finally, for small e, Rse (t) — t- 5 (r) . Thus, since e <C 1, the long-time-scale 
behavior of the time correlation fmiction can be stated as 

< Si (f ) Si (0) >,= 2a,Vee 6 (t) (20) 

Thus, from this point on, we focus on time scales that are on the order of the rocket-dynamics time 
scale, tr, drop the tilde on t, and specify 

< s y, t') s {if, t) >s= 2alTce 5 {ip' - ip) 5 {t' - t) (21) 

as the long-time-scale space-time correlation for s{(p,t), where again we assume delta correlated 
statistics in the angular direction 

Given this correlation, the long-time-scale time correlations for the side load components can 
be finally calculated by combining and averaging instantaneous, non-averaged versions ofEqs. ((§ 
and 

< Fsy{t)Fsy{t) >s=< Fsz{t)Fs, (t) >,= R{t') R{t) AP{t) AP{t) 27Ta^,T^e6{t' - t) (22) 
Obtaining the long-time-scale side load time correlation function represents a crucial step since 
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it allows calculation of side load effects on rocket rotational and translational dynamics. 



5. Pitch and yaw rate response to random side loading 

Given < Fsy{t')Fsy{t) >=< Fsz{t')Fsz{t) >, derived above, we now show that during the period 
when side loads act, < t < T, the long-time-scale evolution of the rocket's pitch and yaw rates 
correspond to Ornstein-Uhlenbeck processes. Given this key result, altitude-dependent means and 
variances for the rocket's pitch and yaw rates, both during the side load period and after, can be 
determined. 

5.1. Derivation of Ornstein-Uhlenbeck rotational dynamics 



As a necessary point of reference, we first note the equations governing rotational dynamics 11 



Im+ {Ixx + l\M\Rl]io^ = (23) 



ICjy + (4^ - /)6<j^a;^ + iuy + \M\ [{L - bf + 0.257?^]^;^ - L_I_^^^t<;^ = Me^t.j (24) 



2 , , \M\LR 



ICoz - (Ixx - I)io^i^y + li^z + \M\ \{L - by + 0.25Ri]u!z + \' j^o ^x^y = ^ext- k (25) 

Regarding notation, moments of inertia are evaluated with respect to rocket- fixed coordinates, 
/ = I{t) and Ixx = Ixx{t) are, respectively, the moments of inertia with respect to either lateral 
coordinate, z or y, and the longitudinal rocket axis, Rr, L, b, and R^, are, respectively, rocket 
radius, rocket length and half-length, and nozzle exit radius, Ve is the speed of the exiting flow from 
the nozzle, M is the associated mass flux, and cOx, ujy, and Uz are, respectively, the roll, yaw, and 



21 



pitch rates. Unit vectors i,j,k, are those of the rocket-fixed system. See 11 for a full description 
and Fig. 4 for a schematic of the rotational and translational variables describing rocket motion. 




Figure 4: Rocket model. Roll, pitch, and yaw rates are shown, respectively, as Ux, ujz, and ujy; 
center of mass position is given by Roit) = [Xo{t) ,Yo{t) , Zo{t)], while rocket-fixed coordinates are 
[x,y,z]. See |11| for further details. 



Since the line of action of the aerodynamic load is assumed coincident with the rocket's longi- 
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tudinal axis, the only force that contributes to the external moment, Me^t, in Eqs. (24) and (25) 
is the in-nozzle stochastic side load. The moment due to these stochastic side loads can thus be 
computed as 

Mext = -{L - b + Xs{t))ix{Fsyj + Fszk) (26) 

As an important aside, and as a means of isolating side load effects, the present model, as well 
as that in does not incorporate wind loading. A simple scaling argument (see Appendix B) 
indicates that random winds: i) under most conditions, do not excite rotational motion, and ii) 
simply function as an additive source of variance in the rocket's translational motion. In other 
words, wind appears to have minimal influence on the stochastic, altitude-dependent evolution of 
rotational dynamics. Rather, (launch-site-specific) mean and random winds simply produce whole- 
rocket, random, lateral translational motion, superposed on a deterministic translational drift. 

During the short side load period, and as found by scaling, it is observed from the equation 



governing roll, Eq. (23), that the coefficients on lo^ and oox, Ixx and (-^^ + 5 M Rj), respectively, 
remain essentially constant. [As shown in |11|, for the model considered, the side load period, T, is 
approximately 10.85 s.] Thus, for an assumed initial roll rate, uJx (0) = 0, and during the side load 
period, < t < T, the rate of roll remains zero, ujx= 0. Indeed, as shown in |11|, and due to the 
absence of roll-inducing torques, this feature holds throughout any given simulated flight. 



Thus, we focus on the evolution of yaw and pitch, as determined by Eqs. (24) and (25), respec- 
tively, and observe that these can be placed in the forms: 

duy = -A(t)ujy{t)dt + B{t)Fsz (t) dt (27) 
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du}^ = -A{t)u}^{t)dt - B{t)Fsy it) dt (28) 

where 

A{t) = [i+ \M\{{L - hf + 0.25i?^] / I (29) 

is a positive damping coefficient (positive as determined by a straightforward order of magnitude 
analysis) and 

B{t) = {L-b)/ m 



Next, from, e.g., [39], we note that gaussian, zero mean random functions that are delta correlated 
in time, such as Fsz{t) and Fsy {t) , can be related to Weiner processes via 

B (t) Fs, {t) dt = /D(t) dW, (t) = ^/D(t}[W, {t + dt) - W, it)] (30) 

where Wz (t) is a Wiener process associated with Fsz (t), D(t) is a diffusion coefficient, and where 
a similar expression holds for Fgy (t). Here, using 

r rt+At /■ rt+At 

J J B{r)B{t') < Fsz {f')Fsz (t'j >s df'dt' = jj^ VD{f') ^/D(¥)dWz {t'')dWz{t') = D{t)At 

(31) 

it is readily shown that 



D (t) = {L- hfl ^ (t) R^{t) A^P{t) 27rcTfrce (32) 



where, consistent with dimensional requirements, D{t) has units of t 
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Thus, we arrive at the important result that during the period < t <T, the random rotational 
dynamics of the rocket correspond to Ornstein-Uhlenbeck processes: 



duy = -A{t)ujy{t)dt + ^/D{t)dWz (t) 



(33) 



du}^ = -A{t)ujz{t)dt - ^jD{t)dWy (t) 



(34) 



We discuss in section 7 the physical features underlying the competition between damping, as 



embodied by the terms involving A{t) in (33) and (34), and amplification of pitch and yaw rates, 
as captured by the D{t) terms. 



5.2. Pitch and yaw rate response 



While in many problems, the damping and diffusion coefficients are constant, here they are not. 



Nevertheless, Eqs. (33) and (34) can be integrated (for individual realizations) via an integrating 
factor: 



uja (t) = exp( — I A [s ] ds 



1^ l^Z) (.s')exp?(^ A(s")ds" )| dW [s 



(35) 



where the subscript a denotes either y or z. Given Eq. (35), the mean yaw and pitch rates can be 
determined: 



< Ua (t) >= exp (- ^ |^I)(s')exp(^ A(s")ds" )| < dVF (s') > 
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or since < dW ( s ] >= 



< UJa (t) >= 



0<t<T 



(36) 



[Note, we are abusing standard notation by expressing the expectation over Wiener processes, 
Eu},s'W {s' + ds') — W{s') as < dW{s') > .] Likewise, (35) allows determination of the time- 
dependent variances for pitch and yaw rates: 



< oolit) >= exp (-2 A{s')ds') /g* {y/ Dis')D{p')exp (// A(.")ds" ) exp [j^' Aip^dp^^ ) } < d/3, 

f^f^{F{s',p')}<dW{s')dWip')> 



s'p' > 



where < dPs'p' >=< dW {s') dW {p') > . Expressing the term /q /J {F{s\p')} < dW {s') dW {p') > 
as a discrete double sum and using AWiAWj = At6ij, it is readily shown that 

<a;2 (t) >=exp(^-2^ A(s')ds'^ ^ |d(s') exp(2j^ ^(s")ds" )| ds' 0<t<r (37) 

It is also easily shown that in the case of constant A and D, this formula leads to the well-known |39| 
variance expression for constant coefficient OU processes. 

During the post-side load period, T < t < Tf, the (short-time-scale) side load terms in equations 



(27) and (28) no longer appear, and thus neither do associated (long-time-scale) Wiener processes 



in Eqs. (33) and (34). Thus, carrying out a single realization integration of the latter two equations. 
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we obtain 



(t) = tOa (T) exp 



Aiq)dq 



(38) 



where coa (T) is a random initial condition on the post-side load pitch and yaw rate evolution 
equations. Since < uJa (T) > = 0, it is clear that 



< Wq, (t) >= 



T <t<Tf 



(39) 



consistent with numerical results in [11] . 
Likewise, forming < (t) > and using 



< col (.T) >= exp (^-2 A {s') ds'^ |ds' exp(2 j'^ yl(s")ds" )| ds' 



(40) 



obtained from Eq. ( 37 ) , the time-dependent variance of yaw and pitch rates during the post-side- 



load period follow as: 



< oJi (t) >= < ul (T) > exp 



-2£A{q)dq 



T <t<Tf 



(41) 



In closing this subsection, we again note that the rocket's rate of roll, cux, is uncoupled from 



moments produced by random side loads (see Eq. (23)); this is reflected in the solution of Eq. (23): 



{t) = (0) exp[ / F {t') dt'] 
Jo 



(42) 
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where F [t') = [Ixx + f l-^l-^el-^xx • Thus, for an assumed initial roll rate of zero, roll rate remains 



zero, ujx (t) = 0, throughout any given flight. Again, computed results in [Tl] are consistent with 
this observation. 



5.3. Pitch and yaw angle response 



Here, time-dependent pitch and yaw angle variances are determined, as well as the time correlation 
function, {tpci{t')ipa{t)) i where ipa represents either the pitch or yaw angle. Determining the evo- 
lution of pitch/yaw variances provides a further consistency check between theory and numerical 
experiments, and, in addition, provides important physical insight into the large, almost explosive 
growth in translational velocity and displacement variances described in section [6] below. The time 
correlation function is needed in order to determine these lateral translational responses. 

During the side load period, single realization evolution of pitch and yaw angles follows by 



integration of ( 35 ) : 



where 



and 



i^a{t) = / expf{q)Ga{q)dq 
Jo 



f{q) = - / 
Jo 



r' 

y/D{s') exp / 
Jo 



(43) 



(44) 



(45) 



Taking the expectation over the Weiner process, dWa{s'), we find, consistent with numerical 
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experiments, that mean pitch and yaw angles remain fixed at zero throughout the side load period: 



=0 0<t<T 



(46) 



Evolution of pitch/yaw angle variance during < i < T is given by 



t rt 



[rM= / exp/(g')exp/(g)(G„(g')Ga(g)><ig'd(? 
Jo Jo 



(47) 



or, usmg 



ti r-ti \ i-min{ti,t2) 

F{s')dW{s')- F{s'')dW{s'')) = F\s')ds', 
Jo I Jo 



t rt 



JO 



{rM= / / exp/((zOexp/(g) 



/ f-min{q,q') 








exp2 [ 2A{s")ds" 


ds'^ 




Jo 





ds' dqdq (48) 



which, by symmetry, becomes: 











exp 2 [ 2A{s")ds" 


(is'^ dq 




Jo 





(49) 



The detailed form of (49) appropriate to the present rocket model is given in Appendix A. 



Following the side load period, pitch/yaw evolution follows by integration of (38): 



,{t)=cOa{T) f h{s')ds' + MT) 
Jt 



(50) 



where 



/i(s') = exp ( - ^ A(s")ds" 



(51) 
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Taking the expectation again shows that average pitch/yaw angles remain at zero throughout the 
post-side-load period: 

iMt)) = T<t<Tf (52) 



where, by (35), {uJa{T)) = 0. Again, this is consistent with results of numerical experiments in [II]. 



Post-side- load variance follows directly from ( 50 ) 



(i^lit)) = (uUT)) (/^ h{s')ds') + 2 MT)MT)) [Ir Hs')ds') + 

+ (V2(r)> T<t<Tj (53) 



where the correlation {uJa{T)tpct{T)) is non-zero and is given by 



(^^a(T)V^«(T)) = / exp/((7)exp/(r)(G,((7)G«(r))dg 

JO 



(54) 



or 



(a;„(r)Va(r)) = r r r exp f{q) exp /(r)/D(sO ^^(s" )9{s)g{s'') {dWa,{s')dWa{s'')) dq 
Jo Jo Jo 



(55) 



or finally by 



(a;, 



,(r)V'a(r))= r exp /(g) exp /(T) P D{s')g\s')ds'dq T<t<Tj 
Jo Jo 



(56) 



Here, g{s') = exp I L A{p)dp\ , and / and Ga are given, respectively, by (44) and (45). The 
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detailed problem-specific form of (56) is again given in Appendix A. 



6. Asymptotic rocket response to small- amplitude stochastic pitch 
and yaw 

Under conditions where the rocket experiences zero-mean stochastic side loading, we anticipate that 
the Euler angles 4>{t), 0{t), and ip{t) all undergo small, random variations about zero. Under these 
conditions, a leading order asymptotic model of rocket lateral velocity and displacement response 
can be used to validate and interpret the numerically estimated statistics obtained via Model I. 

Thus, considering the full translational equations of motion given in ^llj, we observe that under 
conditions where 




(57) 



for all t, where eo ^ 1, these equations assume the forms: 



MXo = (Pe - Pa)Ae + \M\ve " 0.5Cz) (X^ + + Z^) - Mg + 0{eo) 



(58) 



MYo = Ft^ + F,y + 2\M\{L - b)Lo, + 0{eo) 



(59) 



MZo = Ft9 + Fs, - 2\M\{L - b)ujy + 0{eo) 



(60) 



where 



Ft = {Pe - Pa)Ae + \M\v, 



ex 
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is the total thrust force, and where (/'(t), and ip{t) correspond, respectively, to roll, pitch, 

and yaw angles. Other new terms include the rocket mass, M = M{t), the nozzle exit pres- 
sure and exit area, Pg and Ag, the (mach number dependent) drag coefficient, Cd, the rocket 
cross-sectional area, Ar, the altitude-dependent ambient density, pe, and the rocket center-of-mass 



position [Xo{t),Yo{t), Zo{t)]. See Fig. 4, and refer to 11 for a detailed description 



An important note: scaling, as well as results below show that the thrust terms (involving Ft) 



in the lateral equations of motion, (59) and (60), must be included; due to the magnitude of Ft, 



even small, order e pitch and yaws can produce non- negligible, and indeed, dominant lateral forces. 
In addition, due to the restriction of (j){t) to small magnitudes, the following development does not 
hold when the rocket is given a non-zero roll. Note too another abuse of notation in the use of 
dimensional terms in the order symbols. 

Due to the importance of minimizing random lateral velocities and displacements, and due to 
the miniscule effects of side loads on vertical motion [11] , we focus on rocket dynamics in the lateral 
(i.e., y— and z—) directions. The above equations show that during the side load period, < t < T, 
leading order lateral translational dynamics are determined by the summed effects of three random 
forces: 

a) the nozzle side load, F^^, 

b) the lateral thrust component, zti'VV'aiti produced by random variations in pitch and yaw 
angles, and 

c) a second reaction term, ib2|M|(L — b)u}a±, produced by incremental (random) changes in the 
direction of the mass flux vector, M. 
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[Here, 'ipa+ = the yaw angle, ipa- = the pitch angle, uja+ = w^, the pitch rate, and uja- = ujy, 
the yaw rate.] Following the side load period, only the latter two random forcing terms continue 



to act. 



Letting r]{t) represent either Yo{t) or Zo{t), and focusing in turn on the side load and post-side- 



load periods, one can integrate equations (59) and (60) to obtain single realization solutions. Thus, 
over the side load period. 



ft TP ( n ft M{s') ft 



M{s') 0<t<T(61) 



while over the post-side-load period 



7)(t)=77(r)±2(L-6)^ 



Mis'] 



')ds'± [ 
Jt 



j^ds' T<t<Tf (62) 



[Use Eqs. (59) and (60) to choose the signs in (61) and (62). 



In order to determine single realization displacements, first define 



M (s'] 



ft IVl ) 

Hr, {t; to) = 2{L-b)J^ {s') ds' 

F{t;to) -- 



t. Mis') 



(63) 

(64) 
(65) 
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and integrate (61) and (62) to obtain: 



(t) = [ Gr, {s') ds' ± [ (s; 0) ds' ± [ F{s'; 0)ds' 0<t<T 
Jo Jo Jo 



(66) 



r^{t) = i]{T){t-T)± I H^[s'-T)ds' ± I F{s';T)ds' T < t < Tf (67) 
Jt Jt 



Given Eqs. (61) and (62), time-dependent ensemble average lateral velocities are readily calcu- 
lated: 



< F fsM > i*' ) r 



<v{t)>= I 

/o M{s' 



Mis'] 



^mSf^ds' 0<t<T 
M{s') 



ft M{s') f 
< ^ [t) >=< f, (T) >±2{L- b) -j^^< ^a± {s') > ds'± 



,T W^) T<t<Tf 



From Eqs. I\8h, (36) and (39), it is clear that ensemble average lateral velocity components remain 



zero throughout the entire flight period: 



< 7) (t) > = 0<t<Tf 



(68) 



a result that is again consistent with numerical experiments jlll. 
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Variances of the lateral velocity components likewise follow from Eqs. (61) and (62): 



< -n" {t) >-- 



+ /d /o WW? < > ds'dp' + 

+2(L - 6) /* (Va±(p')^a±(.')) ds'dp' 



(69) 



for 0< t < T, and where it is assumed that pitch and yaw rates are uncorrelated with side loads, 

< UJa± (s') Fsa {p') > = 0. 

Similarly, over the post-side-load period, T <t <Tf, 

<if{t)>= < ^2 (T) >+4{L-bf ^ ^ < a;„± (s') a;„± (^ ) > ds'dp' + 

+ /t /t < V'a±(s')V'a±(p') > ds'dy + 

+2(L - 6) W) (V;,i(y)a;„±(s')) d.'dp' (70) 

where we assume that < i) (T) ^^\l^j i^a± (s') ds' > = 0. See Appendix A for problem-specific 



versions of ( 69 ) and ( 70 ) . 



Ensemble averages and variances for lateral displacements can be obtained using the single 



realization solutions in Eqs. (66) and (67). Thus, it is again readily shown that average lateral 



displacements are zero both during and after the side-loading period: 



<n{t)>= [ < Gr, {s') >ds' ± [ <H^ {s- 0) >ds' ± [ F{s')ds' = 0<t<T (71) 
Jo Jo Jo 
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<7]{t) >=<'n{T) > {t-T)± [ < Hr,{s';T) >ds' ± [ F{s')ds' = T<t<Tf (72) 

Jt Jt 



which is again consistent with ^llj. Note, G^, Hrj, and F are given by Eqs. (63), (64) and (65), 
respectively. 

Variances are hkewise found: 



< r/2 (t) >= /J < is') G^{p') >ds'dp' + < Hr, {s'- 0) Hr,{p'; 0) >ds'dp' + 



+ /d /o < F{s'- 0)F{p'; 0) > ds'dp' 0<t<T 



(73) 



< ri^ (t) >= < 7)2 (T) >{t- Tf + JtJt< is'; T) Hr, {p'; T) >ds'dp' + 

+ f^f^< F{s'; T)F{p'; T) > ds'dp' T <t<Tf (74) 

where it is assumed that all cross-correlations between G^, Hr,, F, and i](T) are zero. Problem- 



specific versions of (73) and (74) are given in Appendix A. 



7. Results and discussion 

As an important preliminary, we list essential explanatory remarks and observations. 

A) As shown in Appendix A, most working variance formulae involve one or more time-dependent 
terms: AP{t), R{t), M{t), and I{t). Here, these (non-stochastic) terms are determined using 
numerical data from execution of Model I |11|. In cases where the current model (Model 
II) is used to predict the dynamics of actual rockets, R{t), the time-dependent nozzle radius 
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corresponding to the instantaneous mean separation line position, Xs{t), represents the most 
difficult-to-determine parameter. At least two straightforward approaches, i.e., methods not 
requiring high-level modeling and computation, are available. The first combines one 
of several semi-empirical separation pressure correlations, see, e.g., flS^, with a model of 
nozzle flow upstream of separation. See Fig. 9 in [12]. The second uses scaling to obtain 
an approximate model of mean separation line motion as a function of time (altitude). This 
approach will be described in a future paper. 

B) The parameter AP{t) = Pi{t) — Pa{t), requires estimation of Pi{t), the nozzle- wall pressure 
near the incipient separation point. Again, an alternative to complex modeling and numerics 
can be found in a semi-empirical approach outlined in ^12j . 

C) Examination of the working formulae used to compute theoretical variances (see Appendix A) 
shows that the three a priori unknown boundary layer parameters, (Jg, e, and Tc, representing, 
respectively, the nominal length of the boundary layer-shock interaction zone (fx^ w L^), the 
ratio of the boundary layer to rocket dynamics time scales, and the correlation time for local 
separation line displacements, only appear as the product, Kq = (r'^Tce. Thus, Model II has 
only one available adjustable parameter, Kq, for fitting theoretical to experimental data. 

D) Rather than employing any of a number of standard fitting procedures for estimating Kq, we 
use the following simple approach. First, express Ko as Kq = PoC^o^oTco, where Po becomes the 
fitting parameter and the last three parameters are assigned nominal, empirically- or scaling- 



based values: ao = 2.54 (10 ) m 13 31 , e = 0.1 (scaling), and Tc = 10 s 31,32 . Second 



use straightforward trial and error to estimate /3o- Since the /3o thus obtained is nearly equal 
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to 1/4, we arbitrarily set j3o = 0.25. 

E) It is important to note that if we reasonably (though arbitrarily) take Tc as the most uncertain 
of the three parameters cr, e, and Tc, i.e., express Tc as Tc = Po'^co-, and interpret Tc as the 



approximate frequency of large-scale separation hne motion 34 , then the estimated value of 



^ = 400 Hz, is comparable to frequencies (~ 300 Hz) experimentally observed in shock- 



separated flat plate boundary layers 31 . 



F) As noted in Appendix A, scaling, as well as numerical results (Model I), show that during the 
side load period, < t < T, variance growth in lateral rocket displacements, Yo{t) and Zo{t), 
remains negligible relative to that observed during the post-side- load period, T < t <Tf. For 
simplicity, during < t < T, theoretical (Model II) displacement variances are thus set to 0. 

G) Parameters defining the model rocket can be found in fll|. 

H) An ensemble of 100 numerically simulated rocket ascents are used to estimate experimental 
(Model I) means and variances. Since time-dependent statistics estimated using a smaller set 
of 40 ascents differed by no nore than 14 % from these, no attempt has been made to examine 
larger ensembles. 

I) As noted in section 6, theoretical time evolution of the following mean values: pitch and yaw 
rate, < ujzit) > and < ujy{t) >, pitch and yaw displacement, < 6{t) > and < tlj{t) >, lateral 
velocity, <Yo> and < Zq >, and lateral displacement, < Yo{t) > and < Zo{t) >, all remain 
identically zero throughout the simulated flight period, < t < Tj. These results are, in every 
instance, consistent with the numerical observations in 111]. 
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In discussing the results, we follow three threads: 

A) Evolution of pitch and yaw rate variances turn out to play a central role in evolution of not only 
pitch and yaw displacement variances, but also in variance evolution of lateral translational 
velocities and displacements. In the latter cases, thrust components produced by random 
pitch and yaw emerge as the dominant mechanism generating translational stochasticity. 

B) We find that physical interpretation of many of the results below can be usefully framed in 



terms of the OU damping and diffusion coefficients, A{t) and D{t), given respectively by (29) 



and (32). The utility of these parameters traces to the central role played by random pitch 



and yaw rates in stochastic rocket dynamics. 

C) The transparent physical content of A{t) and D{t) allows straightforward identification of 
practical approaches for minimizing or mitigating against the effects random pitch and yaw. 

7.1. Pitch and yaw rate response 

Numerical and theoretical pitch and yaw rate variances, < uj'^{t) > and < 0Jy{t) >, are compared in 
Fig. 5. General theoretical expressions, applicable over the side load period and post load period. 



are given respectively by Eqs. (37) and (41); corresponding expressions specific to the present 



model are given by (A-1) and (A-3). It is clear that theoretical pitch and yaw rate variances 



remain qualitatively consistent with numerical estimates throughout the entire flight period. 
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0.08 




time (s) 



Figure 5: Comparison of yaw, Wy, and pitch, cj^, rate variance evolution as estimated using Model 
I (solid line) and as computed via Model II (dashed line). The side load period ends at t = 10.85 s. 
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7.1.1. Pitch and yaw rate response: side load period 

During the side load period, theoretical and numerically observed stochasticity, indicated by the 
variances, clearly increase in a qualitatively consistent manner. [Note, differences between ensemble- 
averaged pitch and yaw rate variances decrease when the ensemble is enlarged from 40 to 100 
simulated flights. We anticipate that the observed differences in figure 5 become negligible, as they 
should by symmetry, as — t- oo. A similar comment applies to all other numerical results.] 

Physical insight into variance growth during the side load period, and identification of design 
approaches for minimizing growth of stochastic pitch and yaw, follows by looking at a simpler 



model. Thus, take the damping and diffusion coefficients, A{t) and D{t), in (29) and (32), as 
constant. Although an approximation here, this simplification becomes increasingly accurate as 
the side load period shortens. Under this approximation, variance growth follows the standard OU 
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< ^Ut) >= ^ (1 - (-2^ot)) (75) 



with the associated rate of variance growth, valid in the limit AqT << 1, given by 



dt 



Do + 0{AoT) (76) 



Here, Ao = A{<d) and Do = D{0). 



We use (75) to examine physical features underlying pitch/yaw rate variance in the limit. 



AoT >> 1, appropriate under conditions where damping is strong and/or ascent is slow. Equation 



(76) is used to examine the opposite limit, AoT << 1, i.e., the limit applicable to the present 
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numerical [11 1 and theoretical models. 

Considering first AqT >> 1, the asymptotic variance takes the form: 



' 2Ao /(0)|M| ^ ^ 



where, as determined by scaling, the relationships, |M|(i — 6)^ >> / and |M|(L — b)^ >> 0.25-Rg, 



have been used; refer to (29) 



Physically, the right side of (77) represents the ratio of stochasticity amplification via side-load- 
induced random torques, to stochasticity stabilization via mass flux damping. The former feature 
becomes apparent when we restore the squared moment arm, {L — b)^, cancelled in obtaining (77). 



The role of mass flux damping can be ascertained by sketching the incremental (vector) change, 
AM, in the mass flux vector, M, produced by an incremental pitching or yawing displacement 
(about the pitch or yaw axis). Such a sketch shows that AM acts (nominally) in the plane of the 



nozzle exit, in a direction opposing the angular motion. Thus, with reference to (29) , (33), and 



(34), we recognize that AM produces a retarding, i.e., stabilizing torque that opposes all pitch 



and yaw motions. 



Based on this interpretation, we can use (77) to identify practical approaches, appropriate when 



AqT >> 1, for either reducing random pitch/yaw amplification, or increasing mass flux damping: 

a) reduce the nozzle size (characterized by -R(O)), 

b) reduce the near-exit, ground-level pressure difference, AP(0), e.g., by reducing the degree of 
over-expansion. 
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c) increase the rocket's moment of inertia, as characterized by /(O), and/or 

d) increase the mass flux, m. 

These observations presume, reasonably, that the boundary layer separation parameters, as, Tc, and 
e remain, in an order of magnitude sense, independent of nozzle size and in-nozzle flow conditions. 

Turning to the limit where the side load period is too short for the asymptotic regime to set in, 
AoT << 1, again the limit appropriate to the present study and that in fTTl, we expose features 



that determine the rate of variance growth. Again using initial values for the parameters in (32) 



Eq. ([Tej) yields (to 0{AoT)) 

d < bjlit) > 



Do 



dt 

= (L - 6)2/-2(o)i?2(o)A2p(0)27rcj2Tce 



(78) 



Again assuming that the parameters cJs, Tc, and e remain nominally independent of nozzle size 
and in-flight nozzle flow conditions, we observe that the rate of pitch/yaw rate variance growth can 
be minimized, e.g., by: 

a) again using high-rotational inertia rocket designs, 

b) moving the center of mass toward the nozzle exit (thus reducing the moment arm, L — b), 
while maintaining high I, 

c) using smaller (radius) nozzles, and of course when feasible. 
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d) operating in an underexpanded condition, 

where d) eliminates separation and thus side loading. 

Other approaches, designed to reduce o"s, Tc, and/or e, might include boundary layer manipula- 
tion via, e.g., wall fluid injection or suction, and/or active or passive mechanical forcing. 

7.1.2. Pitch and yacw rate response: post-side-load period 

Side loads cease when the rocket reaches an altitude where ambient pressure drops below the 
near-exit nozzle pressure. From this instant onward, and as is apparent in both numerical and 
theoretical plots of < > in Fig. 5, rocket pitch and yaw, become subject only to the damping 
torque produced by incremental changes in M. 



Equation (41) shows that the post-side-load decay in pitch/yaw rate variance is exponential. 



with the rate of damping increasing with time. Physically, and as shown by (29), accelerated 
rotational damping, in turn, reflects the inverse dependence of the damping coefficient, A{t), on the 
time-decaying moment of inertia, I{t) : as the rocket becomes decreasingly resistant to pitch/yaw 
rotations, the damping torque can effect ever-larger influence on pitch and yaw. 



We can use (29) to quickly identify design strategies that, e.g., enhance post-side-load damping 
of random pitch and yaw. Thus, assuming that the term in is small, the case here, then 

a) increasing M, and/or 

b) increasing the moment arm, L — b, 
will improve damping. 
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Clearly, trade-offs are required since, for example, strategy b) enhances variance growth during 
the side load period. By contrast, increasing nozzle mass flux, at least up to magnitudes for 
which |/| remains smaller than |M|(L — 6)^, is always beneficial since this enhances mass-flux 



damping. [Interestingly, for large enough |M|, the negative term, /, in (29) can become dominant, 
transforming A{t) into an amplification coefficient. A simple analog that explains this effect can 
be found in the gravity-driven pendulum: shortening the pendulum length during motion increases 
the amplitude of the motion.] 

7.2. Pitch/yaw displacement variance 

Pitch and yaw displacement variances are compared in Fig. 6. General theoretical expressions for 



the side load and post side load periods are given respectively by Eqs. (49) and (53), and problem- 
specific versions are derived in Appendix A. As in the case of pitch/yaw rate variances, reasonable 
consistency between theoretical and numerical results is again observed. 

In order to interpret variance growth during the side load period, we again note that due to the 
brevity of the side load period, A{t) and D{t) can be approximated as constant and equal to, say. 



Ao = A{0) and Do = D{0). Starting with (|49|), noting that AqT « 1, and expanding (|49|) in Aot 
then leads to 

<€{t) >- ^t' + o{Ay) 



where 1 >> >> A'^t'^, and where Do is given by (78). Thus, the observed growth in pitch/yaw 



variance during side loading reflects the dominant effect of diffusive, i.e., stochastic, side-load-driven 
growth over mass-flux damping. Variance growth is slow since the diffusion coefficient. Do, is small; 



45 



a quick order of magnitude estimate, yielding Do ^10 ^ s ^, shows that Dot^ /?> is on the order of 
observed and predicted side load period variances. 




0.2r 




Figure 6: Comparison of pitch, 0, and yaw, ip^ displacement variance evolution as estimated using 
Model I (solid line) and as computed via Model II (dashed line). The side load period ends at 
t = 10.85 s. 
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During the post-side-load period, and based on (50), it proves useful to interpret the term, 



fj, h{s')ds' , appearing in (53) as a response function; here, the function yields the total random 
change in pitch/yaw displacement, over the interval from T to t, produced by the random initial 



pitch/yaw rate, oJaiT). Thus, examining the three terms on the right side of (53), we observe that 
the displacement variance at t, < '4)'^{t) >, corresponds to the sum of: i) the average squared 
response to the random input 0Ja{T)., ii) the weighted average linear response to u:a{T), where 
weighting is with respect to the random initial displacement, ipaiT), and iii) the initial displacement 
variance, < V'a(^) > • While the response function increases with time, t — T, it is readily shown 
that the rate of increase decays with increasing t — T; thus, in Fig. 6, we observe a roll-off in 
displacement variance. 

It now becomes apparent that post-side-load displacement variance at any time t increases 
(decreases) with any increase (decrease) in any of the above inputs. This picture provides an 
explanation, for example, of the larger post-side-load numerical (Model I) yaw variances observed in 
Fig. 6. As another example, under circumstances where damping, A{t), decreases, we expect, based 
on simple physical intuition, that the lag between angular velocities and displacements decreases 
and thus, the correlation < LOa{T)ipa(T) >, increases. Thus, consistent with our intuition, weakly 
damped rockets arrive at t = T with a greater range of random pitch/yaw displacements; this initial 
input, in turn, leads, via the response function, to larger post-side- load displacement variance. 
Clearly, practical design considerations, similar to those discussed above, can be identified and 
used to manipulate the response function. 
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7.3. Lateral velocity and displacement variance 

Lateral velocity and displacement variances obtained via Models I and II are compared in Figs. 7 
and 8, respectively. Full theoretical velocity expressions applicable to the side load and post side 



load periods are given respectively by Eqs. (69) and (70). Corresponding expressions for theoretical 



displacement variances are given in (73) and (74). 



The most important observation, as mentioned above and as discussed in Appendix A, concerns 



the post-side-load emergence of the thrust force term: as shown by Eqs. (A-10) and (A-13), post- 
side-load lateral velocity and displacement variances are overwhelmingly determined by the thrust 
force component acting in the random pitch and yaw directions. 

Insight into variance growth over the entire flight period can be gained by interpreting and 



scaling individual terms contributing to the velocity variance in (69) and (70). Thus, taking the side 



load period first, the terms on the right side of (69) correspond, respectively, to velocity variance- 



production via: i) weak random side loads (weak relative to the thrust. Ft)-, ii) translational reaction 
to weak mass flux damping of random pitch and yaw rates, iii) laterally-acting components of the 
thrust force, and iv) coupling of effects ii) and iii). While Ft is large, angular displacements during 
< t <T, relative to those appearing during the post-side-load period, remain small. 

During the post-side-load period, the side load contribution, i), drops out and a contribution, 
< ifir) >, capturing the lumped, integrated effects of ii) through iv), replaces it. During this 
period, angular displacements grow large enough for large lateral thrust force components to appear, 
generating, in turn, large lateral velocities and displacements. 

Notice that lateral displacements, on the order of several kilometers, are observed. As discussed 
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in 11 , scaling readily shows that experimental and theoretical displacements are consistent with 
the magnitudes of computed side loads, thrust forces and angular displacements. Note too, that 
these displacements represent small fractions of the rocket's vertical displacement; the difference 




Figure 7: Comparison of lateral rocket velocity variance evolution as estimated using Model I (solid 
line) and as computed via Model II (dashed line). The side load period ends at t = 10.85 s. 
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between ideal, zero-side-load three-dimensional rocket displacements (relative to the launch loca- 
tion) and displacements observed when side loads are accounted for, remains less than 0.1 % over 
the simulated flight period, (OjTy] [II] . 




Figure 8: Comparison of lateral rocket displacement variance evolution as estimated using Model I 
(solid line) and as computed via Model II (dashed line). The side load period ends at t = 10.85 s. 
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8. Summary and conclusions 

A set of analytical models are proposed which connect stochastic, shock-induced boundary layer sep- 
aration in over-expanded nozzles, to random nozzle side loads and associated rocket response. The 
broad objectives center on establishing the consistency of Models I pTj and II, and on broadening 
understanding of the physical features connecting boundary layer separation and rocket response. 

The physical consistency of Model I is established by demonstrating that the simple model of 
stochastic separation line evolution proposed in [IT] allows direct derivation of the ad hoc side 
load model assumed in [TT], as well as derivation of observed side load amplitude and direction 



probability densities [18 -20 ; the demonstration thus ties the most uncertain elements of Model 
I, the separation line and side load models, to experimental observation. The consistency of the 
analj4;ical models proposed here, collectively Model II, is established by showing that predicted 
altitude-dependent statistics of rocket rotational and translational motion are, in every case exam- 
ined, consistent with those predicted by Model I. 

From a practical standpoint, the analytical models proposed here provide a simple framework 
for analyzing stochastic rocket response to side loading, as well as identifying design strategies that 
minimize either side loads, side load-induced torques, and/or rocket response to these loads and 
torques. 

Appendix A: problem-specific variances 

Detailed variance expressions for pitch and yaw rates and angles, and lateral translational velocities 
and displacements, appropriate to the set of conditions used in numerical experiments, are given 
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here. A few non-trivial derivational details are included. 



Variances during the side load period, < t < T, depend on the diffusion coefficient, D{t), 



and damping coefficient, A{t), given respectively by Eqs. (32) and (29). These in turn depend 
on the time dependent position of the shock in the nozzle, Xs{t) (as given implicitly by R{t)), the 
time-dependent pressure difference, AP(t), the time-dependent rocket mass, M(t), and moment of 
inertia, I{t). Here, these terms are determined via quadratic fits to data from numerical experiments; 
under, e.g., real flight conditions, these terms could, in principal, be estimated using computational 
fluid dynamics simulations, perhaps combined with experimental measurements. 



Pitch and yaw rate variances 



Pitch and yaw rate variances during the side load period, given by Eq. (37), assume the specific 
form: 



-'o 



iHs') 



\ "-^Y'ds' 0<t<T (A-1) 



where 



Cl = 2 



i+\M\{L-bY 



and where, due to the fixed magnitude of M, 



M{t) = Mo- \M\t 



(A-2) 



m 



h - \i\t 
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Derivation of (A-1) is straightforward, requiring separate evaluation of the three integrals in (37) 



with, e.g., the equation immediately above used for the I{t) term in the diffusivity, D{t) (equation 



(32)), and by noting throughout that M and / are constant. 



During the post-side-load period, T < t < Tj, Eq. (41) assumes the specific form: 



Iq - \I\t 
In - \i\T 



Cl 



T <t<Tf 



(A-3) 



Pitch and yaw variances 



The general forms of the side- load-period and post-side- load pitch/yaw variances are given respec- 



tively by (49) and (53). The former takes the specific form 



Iq - \I[ 
lo 



ci/2 



Iq - \I\ 
lo 



ci/2 



dq\ dq 0<t<T (A-4) 



where 



Fi{q) = 2t^ alT,e{L-bf [ 

Jo 



2 fi R^{s')A^P{s' 



lo - \I\S' 



Cl 



ds' 



(A-5) 



Again, this expression follows via straightforward evaluation and integration of individual terms in 



(49). 



Pitch/yaw variance during the post-side-load period is most easily determined by computing 



individual terms in ( 53 ) . Thus 



t N 2 



{co\T)) / h{s')ds 



T 



{lo - |/|r)^2+i - (lo - \i\tr+' (A-6) 
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where C2 = ci/2 = (I +\M\{L - bY) /\I\. In addition, 



* \ 11 



2 {cJa{T)t^a{T)) / h{s')ds' = 2 {uJa{T)ipa{T)) - 

\ ay IH^ay J/\ 1^ y J J \ ay ,^ay ^/ |/|(1 + C2) L 



(Jo - \i\Tr^+^ - (lo - \i\tr+' 

(A-7) 



where the detailed form of {uja{T)ipa{T)) follows from (56): 



(a;,(r)V'a(T)) 



r 


'lo - \i\q 


C2 


lo - \i\T 


7 


Jo 


lo 




h 


Jo 



D{s')g^{s')ds'dq 



and where, from (32) and the definition of g immediately following (56) 



D{s')g\s') = 2TTaiTce{L - b) 



2 , i^2(g/)^2p(^/^ 



Pi 



lo 



lo - \I\S' 



Cl 



(A- 



(A-9) 



The final term in (53 ), < ■0q,(7') >, follows by setting t = T in (A-4). 



Lateral translational velocity and displacement variances 

Full expressions for lateral translational velocity variances during and subsequent to the side load 



period are given respectively by ( 69 ) and ( 70 ) . Corresponding displacement variances are given by 



(73) and (74). Considering the entire flight period, < t < Tj, we find via scaling and numerical 
experiments that in the present model, velocity and displacement variances remain small during 
the side load period, but grow explosively during the post-side-load period. 

Due to the dominance of post-side-load variances, for simplicity, and as mentioned, when com- 
puting theoretical velocity and displacement variances, we simply compute < 'ff'{t) > and < ry^(t) > 
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during the post-side-load period, using only the terms involving the thrust force, Ft- It is readily 
shown that these terms are, at minimum, two orders of magnitude larger than all other, neglected 
terms. 

Thus, the post-side-load velocity variance assumes the form: 



t rt 



Ft{s')Ft{p' 



T JT 



M{s')M(p' 



■JT (V'a±(s')V'a±(p')) ds'dp 



(A-10) 



The time correlation function, in turn, is given as the weighted sum of three separate variances: 



(V-,±(s')Va±(p')> = kis')kip') (ojl^iT)) + {^^^l^{T)) + [fc(s') + Hp')] (a;„±(r)^„±(r)) (A-11) 



where 



1 



1 



iio-\i\TY^\m+c2) 



I,-\I\Tr+^-iIo-\I\q) 



C2 + 1 



(A-12) 



The term (^uj'^_^(T)'j follows by using t = T in ( |A-ip , < 4^1{T) >, follows by doing the same m 



(AJ), and {LOa±{T)tPa±{T)) is given by (jg. 

Finally, lateral displacement variance during the post-side-load period is given by 



t rt 



{v\t))= / / {Fis')Fip'))ds'dp' 

JT JT 



(A-13) 



where 



r FTis')FT{s) 



T JT 



M{s')M{s) 



{'^a±{s)'4}a±{s)) ds'ds 



(A-14) 
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Appendix B: Random wind effects 

When considering the dynamic effects of wind loading, it is useful to focus on the ratio of the 
characteristic time scale associated with say a turbulent cross-wind, = l^/v^, relative to the 
rocket dynamics time scale, tr = Lr/Vr : 



T}^ = h^Y^ (B-i) 

TR Lr 



where and are the integral length and velocity scales for the cross-wind, Lr is the rocket 
length, and Vr is the characteristic rocket speed. Given reasonable, readily estimated velocity scales, 
then wind loading produces non-negligible pitch and yaw responses only under conditions where 
the length scale ratio produces time scale ratios of order one or less. Under these circumstances, 
the turbulent wind loads act fast enough to produce non- negligible pitch and yaw moments. In 
contrast, when the turbulence time scale is much larger than the rocket dynamics time scale, then 
lateral wind turbulence merely produces random lateral displacements, with minimal pitch and 
yaw (where the random translational response presumably has zero mean relative to the mean 
displacements and velocities produced by the mean cross-wind). 

Based, e.g., on a large number of altitude-dependent wind measurements (taken under nominally 



normal daytime conditions at Cape Kennedy) 40 , = 0(10 m/s), while ly^ = 0(10 m), where 



the latter corresponds to the approximate correlation length scale. Thus, since Vr = 0(10'^ m/s) 
and Lr = 10 m, t^/tr = 0(10^ ). While prevailing winds thus have negligible effect on representa- 
tive rocket rotational dynamics, the wind nevertheless produces a random, translational response. 
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Indeed, a quick estimate shows that the lateral wind load (= 0{paV^Aji), where An is the rocket 
lateral area) is of the same order as the characteristic side load (estimated below); however, since 
the side load time scale is fast, specifically, faster than the rocket time scale, tr, then the slow-time 
scale dynamics produced by wind can be neglected when studying the response produced by side 
loads. [When computing the combined effect of random wind and side loading, due to the separa- 
tion in response time scales, the response to the former can be simply superposed on that produced 
by the latter.] 
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